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Abstract 

The entanglement entropy of a subsystem A of a quantum system is ex- 
pressed, in the replica method, through analytic continuation with respect 
to n of the trace of the n th power of the reduced density matrix trp^. We 
study the analytic properties of this quantity as a function of n in some 
quantum critical Ising-like models in 1+1 and 2+1 dimensions. Although 
we find no true singularities for n > 0, there is a threshold value of n close 
to 2 which separates two very different 'phases'. The region with larger n is 
characterised by rapidly convergent Taylor expansions and is very smooth. 
The region with smaller n has a very rich and varied structure in the com- 
plex n plane and is characterised by Taylor coefficients which instead of 
being monotone decreasing, have a maximum growing with the size of the 
subsystem. Finite truncations of the Taylor expansion in this region lead to 
increasingly poor approximations of tr p\. The computation of the entan- 
glement entropy from the knowledge of tr p\ for positive integer n becomes 
extremely difficult particularly in spatial dimensions larger than one, where 
one cannot use conformal field theory as a guidance in the extrapolations to 
n=l. 



1 Introduction 



The entanglement entropy has become a privileged measure of bipartite entan- 
glement for pure states, that is of the entanglement between a subregion of the 
system and the rest of it. In order to compute it, one divides the system in the 
state l^), into two complementary subsystems A and B. The reduced density 
matrix of the region A is obtained by tracing out the degrees of freedom inside 
the region B, pA = tr #|\l/)(\E r |- The correlations between A and B give rise to a 
mixed state for pa whose von Neumann entropy Sa = — tr pa log pa is the entan- 
glement entropy of the region A (see e.g. pQ for reviews). Other quantities such 
as the Renyi (R(a)) [2] and the Tsallis entropies (T(a)) [3] 

R(a) = l0 f tlp % ; T(a) = tTpa - 1 , a>0, (1.1) 
1 — a 1 — a 

are strictly related to the entanglement entropy of the block A. Indeed both of 
them coincide with it in the limit a — ► 1. 

Studying directly the entanglement entropy in the ground state \^f) of a chosen 
Hamiltonian requires the ability to obtain such ground state. The complexity of 
this, in general, grows exponentially with the size of the system (recently new 
methods have been proposed to approximate the ground state of two dimensional 
systems that only scale polynomially with the system size [HE]). 

However, when a in Eq. (II. 1|) is a positive integer n, the Tsallis and Renyi 
entropies can be computed in the context of quantum field theory (QFT), without 
the explicit knowledge of the ground state. The procedure is known as the "replica 
method" [6j [7] and is based on the fact that the expression trp A can be written 
as a partition function of the model on a n— sheeted Riemann surface (or its 
generalization in d + 1 dimensions), divided by Z n , where Z is the partition 
function of the original system . The entanglement entropy is obtained from 
it through an analytical continuation from positive integer n to real a using 
Sa = — lim a _*i ^trp^. Even if all values of tr//| were provided for n G N, 
analytic continuation to real a might not be uniquely determined, as examples 
on disordered systems show [9 J . It is therefore desirable to have some information 
on the analytical properties of tr p°^. 

In the case of 1 + 1 conformal field theories (CFT) the behaviour of the T(n) 
when A is composed by a single interval of length t and is a part of an infinite 
gapless system is known to be [H [7] 



(a\ c(n-l/n)/6 
„ (jj • (1.2) 



Here a is an ultraviolet (UV) cutoff, c is the central charge of the associated CFT 
and r n is a non-universal function of n which is known only in some special cases. 
In some cases even the full analytic continuation has been performed [8] allowing 
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one to extract also the non universal additive corrections to the scaling of the 
entanglement entropy . 

In d > 1, however, CFT results are not available and one thus needs to 
carefully study the analytical properties of tr p°^ before trying to perform any 
analytic continuations of it. 

The purpose of this work is to perform such a study by exploring the analytic 
properties of tr function of a in the whole complex plane for a class of 

quantum Ising-like critical models in both 1+1 and 2+1 dimensions. By doing 
this we uncover a highly non-trivial analytical behaviour of tr p^ for a real in the 
interval < a < 2. In all the models considered, we find a threshold value n c in 
the interval 1 < n c < 2 which divides tr//| into two 'phases'. In the region with 
a > n c the quantity tr p^ behaves very smoothly as a function of a; by expanding 
it in a Taylor series centred around any of the a a > n c , we extract alternating 
series with monotone decreasing coefficients, then, any truncation ^fc=o c fc( a ~~ 
a Q ) k can be used to approximate tr p^, the error of this approximation being less 
than \c n (a — a ) n \. For this reason we call the region a > n c 'perturbative' region. 

On the contrary, in the other 'phase', the Taylor expansions of tr//| with 
< a Q < n c has coefficients whose moduli are no longer monotone decreasing. 
They show a characteristic peak whose height increases with the size of A (see 
for instance FigfT]). In this regime any truncation of the Taylor series becomes an 
increasingly poor approximation of tr p\ as the size of A grows. Thus we name 
the region < a < n c 'non-perturbative' region. 

We shall see that the threshold value n c coincides with the point where the 
expression da - 2 A + g da3 A changes sign, that is 



In the models we have considered it turns out that the two different 'phases' 
are separated by a smooth cross-over in both d = 1 and d = 2. However we 
cannot exclude that, in higher dimensions, this cross-over turns into a real phase 
transition. In fact it has been recently observed in Ref. [10J that the Renyi 
entropies in the 0(N) model undergo a phase transition as a function of a for d 
close to 3. This happens in the same range of a in which the models we analyse 
develop the cross-over between the perturbative and the non-perturbative regions. 

The above scenario has important consequences when trying to extract the 
entanglement entropy using the replica trick. On one side the absence of a real 
phase transition makes the analytic continuation (it involves) possible. On the 
other side, the different behaviour of the two 'phases' implies that finding it can 
be really a non trivial task. 

In fact, it is easy to find a simple function f(n) which accurately reproduces 
tr p n A for integer n. However, if we make the 'obvious' extension to real n, by 
assuming that this functional form is also valid for real n, we obtain an incorrect 
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value of the entanglement entropy. 

By exemplifying this we also have the opportunity to compare two different 
methods. One is based on a recent algorithm proposed in Ref. [13] which allows us 
to accurately evaluate the reduced density matrix of a two-dimensional quantum 
system defined on a lattice of small size (in this paper we generally refer to this 
method as the Hamiltonian approach). The other one is based on the replica 
method and applies a simple technique described in Ref. [H] to directly measure 
tr p\ as the vacuum expectation value of a suitable observable in a Monte Carlo 
simulation, in any Euclidean lattice field theory. 

The two methods agree on the Tsallis entropies T(n) for a wide range of 
integer n > 2, once the UV cutoff of the two approaches are matched through 
a suitable conversion factor (this by itself is a non trivial result that provides 
a strong confirmation of the validity of both approaches). However there is a 
mismatch in the results for the entanglement entropy. Its value directly computed 
in the Hamiltonian approach (the correct one) is different from the one obtained 
in the Monte Carlo approach from the naive continuation of the behaviour of the 
Tsallis entropies for integer n to real a. This is a consequence of the intricate 
landscape of trp^ close to n = 1. The scenario we have outlined is common to 
all critical quantum systems in one and two dimensions we have analysed. 

The contents of this paper are organised in several sections. In section [2] we 
describe the quantum models we want to study with different numerical methods 
and, where available, analytical results. In section [3] we describe both the Monte 
Carlo and the variational techniques we use to analyse quantum systems in 2+1 
dimensions. Section O and [5] present the main results of this paper. In Section 0] 
we outline the presence of two 'phases' for tr by considering the features of its 
Taylor expansions. In section [5] we look for the Lee- Yang zeros and other possible 
singularities in the complex plane of a and we find no singularity related to the 
transition between these two regions. We unveil the presence of a smooth cross- 
over between the two regions. In Section[6]we discuss the effects of this scenario on 
entanglement entropy calculations. There we also present the comparison between 
the numerical results obtained with Monte Carlo and variational techniques. We 
conclude with a discussion of the results in Section [7J 



2 The models 

We analyzed the reduced density matrix tr p n A of the ground state near and at 
the quantum critical point in a series of one-dimensional spin \ chain systems 
and their two-dimensional generalizations. In particular, we considered the XY 
Hamiltonian 

H XY = -Y J [(1 + 7K*f + (1 " 7K*j] - 2 h of , (2+) 
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where denote first neighbor sites, erf are Pauli matrices and < 7 < 1, 

h > 0. The subsystem A is a block of neighboring nodes. 

We applied mainly the numerical methods described in [Til E], but we used 
also analytical results. Analytical calculations exist for the critical XX chain 
(7 = 0) [16], for the non-critical XY chain in a field [17] as well as for the 
Heisenberg chain [18] . 

Simplifications arise when mapping the quantum chains into classical two- 
dimensional classical spin systems [19] . The quantum XY chain can be mapped, 
at least for a part of the parameter range, into a classical Ising model on a 
triangular lattice [20] . In this way simple formulas can be easily obtained. For 
instance, if the subsystem A is a block of spins of size L much larger than the 
correlation length £, embedded into an infinite chain with periodic boundary 
conditions at zero temperature in the ordered phase ((i.e. h < 1) the reduced 
matrix is independent of L and can be written in the simple form [21] 



where the parameter e which sets the scale is a known function of 7 and h \20\ [2T] . 
The factor of two in front of (|2.2I) arises from the asymptotic degeneracy of the 
ground state in the broken Z2 symmetry for h < 1. The exponent of two in (j2.2|) 
reflects instead the fact that the segment L in a chain with periodic boundary 
conditions has two points of contact, i.e. two interfaces with the rest of the 
system. 

We are interested in the critical limit e — ► where the infinite products (|2.2|) 
become slowly convergent, therefore we rewrite (|2.2|) in a more suitable form using 
the identity 

TT (1 + (? „ a) = ( _ i7rQT/12 ) T = i L (2.3) 

mar) tt 

71=1 



thus 
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trp? = 2 -9 — / ; v : . (2.4) 
1 77(07-) ry(2r) cv 1 



?7(2aT) 7/(t) 

we introduced the Dedekind rj function 

00 

ri(r)=q^H(l-q n ) (2.5) 



n=l 



because it has a simple behaviour under the modular transformation r — ► — 1/r, 
namely 

= -ri=T7(--) , (2.6) 
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and this leads to rapidly converging products. Precisely we have 



ICi (l + exp(-4n)) Q \ 

n^^i+expc-^)); 

Since ^ ~ log£, the first exponential tells us, according to [HE], that the associ- 
ated CFT in the critical limit is the free fermion model with central charge c=\ 
as expected, however the exact expression (|2.7p allows to discuss the analytical 
structure of the reduced density matrix in the whole complex replica plane. This 
discussion is postponed to Section [SJ 

For the other critical models, to our knowledge there are no analytical results 
available for the reduced density matrices, and we have to rely on a set of different 
numerical techniques. For the generic one dimensional critical XY chains we use 
the expression of the reduced density matrices of an interval in terms of the free 
fermionic modes described in Ref. [11] (for more recent results see e. g. Ref. |12j ) 
and numerically diagonalise it. In the case of finite chains at the critical point 
we use the one dimensional version of the variational technique described in Ref. 
[13J that employs a Tree Tensor Network as an Ansatz for the ground state of the 
system. 

In two dimensions we analyse the Ising model in a transverse field and the XX 
model on finite tori with both the Monte Carlo technique described in Ref. [H] 
and the variational technique of Ref. [13] that uses again a Tree Tensor Network. 
In 2D the subsystem A is one half of a torus bounded by two parallel straight 
lines. 




3 Monte Carlo simulations and variational algorithms 

In QFT the quantity tr p n A can be written as the vacuum expectation value of a 
suitable observable O defined on a larger system composed of n decoupled copies 
of the original system. 

The canonical partition function Z = tr e~@ H of a d— dimensional quantum 
system at inverse temperature (3 can be computed in a standard way by doing the 
Euclidean functional integral in a <i+l-dimensional hyper-cubic lattice A = {x, r} 
(xi,r £ Z) over fields 4>(x) = 4>(x,t) periodic under r — > r + (3. Therefore the 
system composed of n independent replicas of the original system is described by 
the n— th power of Z: 

/n 
[IPWe" E « SH , (3.1) 
k=i 

where <f>k is a field configuration associated to the k— th copy and S[<p] is the 
Euclidean action. 
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In the replica approach to entanglement the subsystem A establishes a process 
of transferring information among the n replicas through a specific coupling: the 
lattice links coming out of nodes of A and directed into the (Euclidean) time 
direction r cyclically connect the copy k with the copy k + 1 [BJ [7] . Let us denote 
with Sa [4>i i 4>2 1 ■ ■ ■ , 4>n] the corresponding coupled action. It is easy to see that 
the quantity 

q _ e ~[S i X ) l'l>i,'p2,...,<t>n]-'E,h=iS[4>h]) ^ 2) 

has the desired property. In fact its vacuum expectation value in the system of n 
independent copies of the original system is 

Z n (A) is the partition function of the system in the n coupled replicas. 

This method can be simply implemented in Monte Carlo simulations. Here 
we apply it to the case in which the system is a periodic lattices of size 1? x 8L, 
where the longer direction is the Euclidean time. We verified that the Euclidean 
time is large enough to ensure the absence of measurable finite temperature ef- 
fects. The Monte Carlo simulations were performed with a particularly efficient 
implementation for the Ising and q— state Potts model described in p3] and with 
the standard action of the isotropic Ising model on a cubic lattice 

S = ~Y J S i S j ; Si = ±l. (3.4) 

(ij) 

All simulations were made at the critical coupling (3 ~ 0.221652. 

With this method we compute directly the Tsallis entropies for integer n. A 
discussion of the results we obtained is postponed to Section [H 

Using the techniques of Ref. [13] we also extract an accurate variational ap- 
proximation to the ground state with a Tree Tensor Network of the corresponding 
quantum model. This is the quantum Ising model on a square torus at its critical 
point obtained by setting the parameters of the Hamiltonian in Eq. 12.11 to 7=0 
and 2h = 3.044 l oj. This approach is limited by the exponential growth of the 
Hilbert space with the system size. The benchmark calculations presented in Ref. 
[13] show that, nevertheless, it can be used safely at least for tori of dimension 
up to 10 x 10. It has the strong advantage that it allows us to extract the com- 
plete spectrum of the reduced density matrix pa for arbitrary blocks. It thus 
provides a way to directly compute tr p a for an arbitrary complex a and also to 
directly compute the entanglement entropy. A detailed explanation on how to 
compute the reduced density matrices in this approach is contained in Ref. |13j . 
We also postpone the discussion of the numerical results for both the Tsallis and 
entanglement entropy extracted with this method to Section [6l 
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4 The two regions of Tsallis entropies 



As we discussed in the introduction the Taylor expansions of tr p^ can be used 
to study the behaviour of the Tsallis and Renyi entropies as a function of a and 
identify two very different regions. In order to exemplify this idea we have to 
express the Taylor expansion of tr p^ about a = 1 in a convenient form. Let pi 
be the non vanishing eigenvalues of pa- At first the normalisation of pa implies 

trp A = ^2pi = l. (4.1) 

i 

Through the chain of identities 

i i k=0 ' k=0 

OO /. \ /. oo 

_ ^ trp A (logp A )' c 1 x fc _v^„ i„ i\k 



( a _l)fc = ^ Cfe(a _l)fc , (4.2) 



^ k\ 

k=0 k=0 

we define the Taylor coefficient c& = tTpA ^°f PA ^ . Since the eigenvalues of tr pa are 
smaller than 1 the Taylor expansion is an alternating series, i.e. Ck+i = — c/~. This 
is an important feature when discussing the possible location of the singularities 
in the complex a plane or in estimating the error we make in truncations of 
the Taylor expansions, as we shall see later. From the spectrum of pa we can 
explicitly compute the above Taylor expansion. 

As an example, the method described in Ref. [Tl] can be used to compute 
the reduced density matrix of a block of spins of an infinite critical Ising chain. 
From it, we compute the moduli of the first Taylor coefficients of tr p^ and plot 
them in Fig{T](a) for four different values of L. It is interesting to notice that the 
plot shows a peak at some value of k > 1 and that its height increases with L. 

A similar behaviour can be observed in a quantum two-dimensional critical 
Ising model on a torus L x L as shown in FigfT] (b). There we plot again the 
absolute value of the Taylor coefficients of the expansion of tr p°^ about a = 1 
for four different values of L. The subsystem A in this case is half a torus. The 
reduced density matrices pa are obtained with the techniques described in Ref. 
|13j . By comparing the same expansion for different systems sizes, L = 4,6, 8, 10 
we appreciate how the peaks move with L and, more importantly, their height 
increases with L. 

By varying the expansion centre of the Taylor expansions a a we can check the 
behaviour of the coefficients in other regions of the a axis. The expression for the 
coefficients of the Taylor expansion about a = a Q is 

r i ^p a A °(\ogpA) k . 
Ck[ao] = ^ • (4.3) 
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Figure 1: Plot of the moduli of the coefficients of the Taylor expansion of of tr p\ 
about a = 1 for (a) the ID Ising model in transverse field at the critical point and 
(b) the 2D Ising model defined on a torus of L x L. In (a) the block A consists 
of L adjacent spins of the infinite chain with L = 64, 128, 256, 512. The spectrum 
of the reduced matrix is obtained applying the method described in Ref . [TT] . 
In (b) the subsystem A is half a torus L x L with L = 4, 6, 8, 10. The spectrum 
is obtained by using the method described in Ref. [13]. In both dimensions the 
height of the peaks of the Taylor coefficients increases with L. 
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We repeat the same study presented in Fig. [I] in the region 1 < a Q < 2 and 
observe that the peaks of Fig. [T] (a) and (b) start moving to left as we increase 
the value of a Q to eventually disappear at a threshold value a = n c < 2. In 
the region where the peaks are present we cannot obtain a good approximation 
of tr by truncating the power series (|4.2p down to only few terms. Also, the 
convergence of the series (I4.2p deteriorates with the size of the system, preventing 
the possibility to find a good truncation scheme independent of the system's size. 
For this reason we call the region a < n c the non-perturbative region 0. 

For a Q > n c the coefficients Cfc[a G ] become monotone decreasing in modulus. 
As a consequence, in this region, the partial sum Sk a = J2k°=o c fc[ a o](« — a o) k can 
be used to approximate trp^ with an error smaller than Cfc 0+ i[a ](a — a ) ko+l . 
In this case thus, we call the region a > n c the perturbative region. 

The threshold value n c is defined as the minimal value of a such that 



that is, the Taylor coefficient Ck are monotonically decreasing functions of k. In 
all the models we have considered, as soon as the criterion in Eq. 14.41 is fulfilled 
by the second and third coefficients, that is, | c$[a]\ < | C2[a]|, it is also fulfilled 
by all the others pairs. Using this property we can define n c as the solution of 
Eq. (|1.3p anticipated in the introduction. We can also use 



the sign of the expression (jl.3p defining Tic i 

sort of order parameter distin- 
guishing between the perturbative (s > 0) and non-perturbative (s < 0) 'phases'. 
As expected, the value of the threshold depends on the size of the system, i.e. 
n c = n c (L). By analysing this dependence we found that in our critical systems 
n c (L) obeys a simple power law 



model. It is also interesting to notice that n Q < 2 for all the models we have 
analysed (see Figs. [2] (a) and (b)). 

5 The nature of the transition 

Having uncovered the presence of two different regions for Tsallis (and Renyi) 
entropies as functions of a, we now want to investigate whether these two regions 

* Starting from here, we will gradually drop the notation a in favour of a. The reader should 
not be confused. Indeed the Taylor coefficient for tr do not depend on a but rather on 
a , the centre of the expansion. However a has the same support than a and we can safely 
simplify the notation. 



Cfc+i[a]| < \ck[a]\ , k = 0, 1, ... . 



(4.4) 




(4.5) 



n c (L) =n a + b/L e , 
where e = i in the ID critical Ising model and e 
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Figure 2: (a) Scaling behaviour of the threshold n c that separates the perturbative 
from the non-perturbative region, as a function of the size L of the block, (a) ID 
quantum critical Ising chain. The rhombi refer to blocks embedded into an infinite 
chain, while the circles refer to half of a finite chains of length 2L. The threshold 
values located respectively at rtj = 1.414(5) and = 1.391(1) are approached in 
both cases with a power law with the same exponent 1/8 . Data for the infinite 
chain are obtained using the method of Ref. [11] while for finite chains with the 
one dimensional version of the method of Ref. [13]. (b) 2D quantum critical 
Ising models on several tori. The blocks are half tori. Here also the location of 
n c approaches its thermodynamic limit as a power law with exponent 1/2 and 
the limit is located at nt = 1.835(1). The data are obtained with the technique 
introduced in Ref. [13]. 10 



are well separated phases divided by a transition point or if instead they are 
different regions of the same phase smoothly connected through a cross-over. Of 
course, these two scenarios have very different physical implications. In the case 
where the two regions are separated by a true phase transition for some real n c in 
the range 1 < n c < 2, there would be no possibility to analytically continue the 
Tsallis entropy from one phase to the other. This, for example, would imply that 
the entanglement entropy could not be obtained via the replica method. This 
is unlikely to happen in one dimensional critical systems since the entanglement 
entropy has been already successfully computed with the replica trick in some 
1+1 CFT [8j. In two dimensions however, for the models we are considering, 
the entanglement entropy has not yet been obtained analytically and one can be 
concerned about if with the replica trick is a viable method to perform such a 
calculation. 

In order to settle the nature of the transition between the perturbative and 
non perturbative regions we start by considering the radius of convergence of the 
Taylor expansion of tr p A around a = 1. It indeed measures the distance between 
the centre of the Taylor expansion and the nearest singularity in the complex a 
plane. As the expansions of tr p A is an alternating series its possible singularity 
belongs to the real axis at the left side of the expansion centre. We can also argue 
on general grounds that the singularity should be at a greater or equal to zero. 

Indeed we expect at least a singularity of tr p a A at a = 0, since tr p A measures 
the number of degrees of freedom of the subsystem A. This number diverges 
exponentially with its size (Eq. (|1.2p . for example, presents an essential singularity 
at a = 0). 

For the sake of completeness one should remember that, when the system is 
in a product state (situation that does not arise close to a critical point) there is 
no singularity at a = 0, even in the thermodynamic limit. This happens, i.e., in 
the XY model described in (12. ip as one approaches the disorder line h 2 + 7 2 = 1 
in the ordered phase (h < 1) . There the ground state becomes the superposition 
of two product states [24] and tr p a = 2 1 ~ a . Similarly, there is a class of z = 2 
quantum critical points in two spatial dimensions where a very simple form for 
tr p A has been proposed [25], namely tr p A = (Zd/Zf) 01-1 , where Zjj and Zp are 
the partition functions of the associated system with suitable boundary conditions 
which select the subsystem A. 

In the cases of highly entangled states we consider here, however, we do expect 
a singularity at a = in the thermodynamic limit and therefore the radius of 
convergence of the expansion (14.2H of tr p A centred around a Q = 1 should not 
exceed 1. A precise value for it can be extracted using the ratio test. This states 
that if the limit r = liim_ 



exists, it coincides with the convergence radius 
r. In a finite system, there is no true singularity and the radius of convergence of 
the series is infinite. However we expect we are still able to identify what would 
turn into a true singularity in the thermodynamic limit, by studying the sequence 
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We anticipate that, notwithstanding there are clues to the presence of singu- 
larities in the thermodynamic limit, a careful analysis allows us to exclude it: the 
only possible singularities of tr p a lie on the imaginary axis of a. 

These ratio sequences show a sort of plateaux at a value r in a certain range 
of k before they eventually run away for k > k Q . The point k Q moves to the right 
within the size of the system, revealing longer and longer plateaux. This has been 
observed in all the models we have analysed. For instance, we may compute the 
required ratios from the Taylor coefficients of a L = 512 block of a critical 
one dimensional Ising model, by using the coefficients plotted in Fig. [I] (a). These 
ratios are plotted in Fig. [3l 

There we can see that the ratios oscillate for a while around a value smaller 
than 1 before running away. This suggests that a new singularity on the right 
side of the one we already expect at a = might exist in the thermodynamic 
limit. 

This scenario is further supported by studying the Taylor expansions around 
other positive values of a. If there were a true singularity on the real axis, the 
function 



would converge to the coordinate location of the singularity in the thermodynamic 
limit and would not dependent on the value of a. In a finite system we expect that 
A[a, k] should approach to a constant value independent of a for a certain range 
of k, and this range should increase with the size of the system. This behaviour 
can be observed, for instance, in Fig. 0] (a) where we plot A[a, k] for the L = 512 
block of the critical quantum Ising chain (the same already considered in Fig. Q] 
(a)), and in FigH] (b) where we plot the same quantity for the two-dimensional 
Ising model on a torus of side L = 6. 

Still all the above signatures are not sufficient to conclude that there is a 
true singularity for real a in the thermodynamic limit in tr p a . As shown by Eq. 
(|3.3p . trp^ can be written as the ratio of two partition functions. This allows us 
to apply a Lee- Yang [26] analysis to it. 




In discrete systems of finite size partition functions are analytic with respect 
to their parameters because they are finite summations of positive terms. For 
the same reason they do not have zeros on the real axis. However there can exist 
zeros on the complex plane of certain control parameter, a in our case, which are 
generally named Lee- Yang zeros [26] . In general, the Lee- Yang zeros are apart 
from the real axis as long as the system size L is finite. However, if a phase 
transition occurs at a critical value n Q £ R, the distribution of zeros becomes 
dense and accumulate on curves. The singularities of the free energy associated 
to this transition lie on the end points of these curves. As L — > oo these end 
points move towards the real axis, eventually touching it at the location of the 




(5.1) 
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Figure 3: (a) Ratio as a function of k obtained from the Taylor coefficients 
of tr/?^ about a Q = 1 plotted in Fig. [TJ The ratios oscillate for k < 20 around 0.8 
before they run away. This could be a hint of a finite radius of convergence for this 
Taylor series in the thermodynamic limit, generated by a conjectural singularity 
located somewhere around 1 — 0.8 = 0.2. In (b) the ratio test shows a similar 
behaviour for the d = 2 critical Ising model on a torus L x L with L = 6. 
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Figure 4: Plot of the function A[a, k] defined in (|5.ip for different values of a. 
This quantity should converge in the thermodynamic limit to the location of the 
alleged singularity. The system in (a) is the same quantum Ising chain of Fig{T] 
The system in (b) is the critical Ising model on a torus L x L with L = 6. Note 
that the position of the flexes on the y axis is nearly independent of a. 
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critical point n . 

In the case of tr p\ = besides the surviving zeros of the numerator there 

could also be poles due to zeros of the denominator. As an example, the exact 
expression of the reduced density matrix in Eq. (12. 1\ for the non-critical XY 
chain, presents a singularity at a. = (the same as the singularity of (|1.2|) ) but 
no zeros. It contains instead a huge number of double poles at a = ±i j^zw e f° r 
any pair of integers k, n £ N. This set is dense on the imaginary a axis. The 
plot of $te tr p2 is presented in Fig. [5] (a) and reflects this intricate structure. A 
similar landscape can be observed by considering the quantum two dimensional 
critical Ising model on a torus (see Fig. [5] (b)) where one can locate few zeros 
(see e. g. Figured]) close to the region Kea = 0.3 where the ratio test give hints 
of a singularity. 

The important outcome of this study is that in all the models we have analysed 
both in d = 1 and in d = 2, these Lee- Yang zeros, which are true singular points 
of the Renyi entropy (jl.lf) in the complex plane, do not become denser as the 
size of the system increases nor they approach the real axis (see Figure [7]) . In 
the thermodynamic limit then, the apparent singularity on the real axis identified 
by the ratio test of Taylor coefficients does not occur. The behaviour of flexes 
described by Eq. (15.ip is probably only a consequence of the presence of some 
structures in the complex a— plane (zeros and/or peaks) which do not evolve into 
a true singularity in the thermodynamic limit. 

The lesson we learn from this analysis is twofold. On the one hand there is no 
true phase transition between the two regions where the Tsallis entropies show 
very different behaviour. The transition is only a smooth cross-over. On the other 
hand the behaviour of tr p a shows an intricate landscape in the region 3f?e a < n c 
while it is very smooth in the complementary region SRea > n c . As a consequence 
a naive analytical continuation of the results obtained for the Tsallis entropies 
in the perturbative region is likely to fail to reproduce the correct entanglement 
entropy. We will provide an example of this fact in the following section. 

As an aside, it is interesting to observe the striking resemblance between 
the location and shape of the zeros distribution of Fig. [7J and those found in 
certain disordered systems [27] . In the latter case the replicas are coupled together 
no longer by the subsystem A, but by the quenched average, which restores 
the translational invariance of the system. As a consequence the latter system 
develops in thermodynamic limit a true singularity associated to replica symmetry 
breaking. 

6 The entanglement entropy mismatch 

In order to test the consequences of the intricate structure of tr p a in the complex 
plane we have performed a numerical experiment. This consists in computing the 
entanglement entropy using the replica trick from the results obtained in a Monte 
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(a) 




(b) 

Figure 5: Plot of Ketr/^ for the non-critical XY chain in the complex replica 
plane (a) and in the 2D Ising model at at the critical point for half a 6 x 6 torus 
(b). Both plots show a very rugged landscape with several dips and peaks. 



16 



trp 




Figure 6: Plot of the real and imaginary parts of tr p a along the imaginary axis 
^sma with Kea = 0.3 for the 2D Ising model on a 6 x 6 torus at the critical 
point. The imaginary part is the curve intersecting the origin. The intersection 
of the two curves identifies a pair of zeros at Qma = ±0.3655. 

Carlo simulation. For the 2D Ising model on a torus indeed there are still no an- 
alytical predictions on value of the parameter appearing in the proposed form 
of the scaling of the entanglement entropy. Some numerical results have been 
recently obtained in Ref. [13]. Here we have computed trp^ for n = 2,3,4,5,6 
with the Monte Carlo technique described in Sect. [3]when A is a half of the torus, 
the same geometry studied in Ref. [13] . The results we obtain agree with the ones 
computed from the data of Ref. [13] provided we introduce a suitable conversion 
factor between the lattice spacings of these two different regularisation schemes. 
A good agreement is indeed obtained by multiplying the lengths measured in the 
Hamiltonian approach by a conversion factor / ~ 1.45 (see FigJH]). This agree- 
ment is a good check of the validity of both methods. To our knowledge it is also 
the first direct comparison between two different non-perturbative techniques on 
entanglement entropy, one involving quantum mechanical tensor network simula- 
tions and the other one using statistical mechanics Monte Carlo integration. This 
is what we expect from universality, indeed, these two methods deal with very 
different systems (in one case the model studied is the 3D classical Ising model 
while in the other case it is the 2D Quantum Ising model). They show the same 
behaviour close to the critical point because they belong to the same universality 
class. 

We can now try to extract a function f(L,n) from the data describing the 
behaviour of tr p n in the whole range of n and L considered . The function should 
depend on only few parameters and is chosen in such a way to agree, in the n — > 1 
limit, with the expected general form of the critical entanglement entropy in 2+1 
dimensions [23j [10] . A sufficiently good fit (a reduced x 2 = 1-75 for 17 degrees of 
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0.3 0.4 0.5 0.6 0.7 

Figure 7: Finite size scaling of the zeros of tr in the complex a-plane. for 
(a) several blocks of different length in the one dimensional XY model and (b) 
several half tori of different size in the 2D quantum Ising model at the critical 
point. In both cases, the zeros do not approach the real axis nor their number 
grows with the size of the system. This fact demonstrates that they cannot 
originate a singularity for real a in thermodynamic limit. 



freedom) is given by the three-parameter function 

f(n, L) = exp[-a (n - -1)2 L - (b + ^-){n --)] 



(6.1) 
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Figure 8: Fit of the Monte Carlo estimations of tr (rhombi) to the function 
(|6.ip in a 3D critical Ising model. The circles are the values calculated in the 
corresponding quantum system in 2+1 dimensions with the method of [13] with 
the length rescaled by a suitable conversion factor. 

with 

a = 0.01177(92) , b = 0.0594(55) , c = -0.378(73) . (6.2) 

The 'obvious' extension to real n, is performed by assuming that this func- 
tional form should also be valid for real n. The entanglement entropy then is 
given by 

7 d 

Sl = — tr pa log pa — - lim — f(n,L) = 8 a L + 2b + c/L. (6.3) 

n-*l an 

Due to our choice of the parametrisation (|6.ip the entanglement entropy has the 
expected functional form Sl = c\ 2L + cq + The numerical coefficients we 
extract however, apart from the 'area term' c%, do not coincide with the ones 
extracted from the direct study of the entanglement entropy computed in Ref. 
[13J. If we call the parameter extracted in Ref. [13] Sl = c / 1 2L + c' + taking 
into account the conversion factor /, we would expect from the values reported 
in (gT2J) that 

c[ = ci / = 0.0682(5) , Cq = c = 0.12(1) , c_ x = c_i// = -0.52(10). (6.4) 

On the other hand their direct numerical determination [13] provides different 
results 

ci = 0.06722(18) , c' = 0.0250(21) , d_ x ~ . (6.5) 
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Figure 9: Plot of the moduli of the coefficients of the Taylor expansion about 
n = 1 of Eq. (jl.2|) with r n = 1. Although it fits perfecltly the behaviour of trp^ 
of the same system of Fig. [1] (a) with L = 512 in the range n > 2, its analytic 
continuation to n ~ 1 region does not reproduce the true Taylor expansion plotted 
in FigH](a). 



We attribute the mismatch of the entanglement entropy to the the different 
behaviour of the the regions 5Rea < n c and $tea > n c . 

A further confirmation of this can be obtained by computing, in the Hamilto- 
nian approach, tr p\ for some real (or complex) a in the non-perturbative region. 
If we add to the fitted data some of these points the quality of the fit to (|6.ip 
rapidly deteriorates. 

One can check that a similar phenomenon holds also for d = 1 models. Indeed, 
considering again the critical XY system in which the subsystem A is composed 
by a block of L = 512 adjacent spins, tr can be evaluated with the method of 
|llj and the data with n > 2 can accurately fitted with Eq. (jl.2p ignoring the n 
dependence of the non-universal coefficient r n which is known to be neglegible in 
such a range |22j . If now we assume that this functional form with r n = 1 is also 
valid for real n and try to evaluate this quantity near n = 1 we obtain incorrect 
results for the entanglement entropy. As an example the Taylor coefficients of 
this fitting function are plotted in Fig. [9] and should be compared with the same 
coefficients directly evaluated on trp^ (see FigfT] (a) ). They show a completely 
different shape. Again this mismatch has to be attributed to having applied to 
the fitting function (jl.2p with r n = 1, which accurately describes the behaviour 
of tr p n in the perturbative region (i.e. n > n c ), the 'obvious' extension to real n 
close to n = 1 which lies in the non-perturbative region. 



20 



7 Discussion 



In this work we analysed the trace of the n— th power of the reduced density 
matrix tr p 1 } associated to a subsystem A in some Ising-like quantum models in 
both one and two spatial dimensions near and at a quantum critical point at zero 
temperature. We unveiled that, depending on the value of n, tr p n A shows two very 
different behaviours. If n is larger than a threshold value n c -located in the interval 
1 < n < 2- its behaviour is very smooth. In particular the alternating coefficients 
of its Taylor expansions, c& are monotone decreasing. In this regime the error 
involved in a truncation of the series does not exceed the first excluded term and 
can be always chosen to be very small. We label this region 'perturbative' region. 

In the range n < n c , on the other hand, the |cfc|'s are no longer monotone 
decreasing function of k but they present a characteristic peak (see Figfj] ) whose 
height increases with the size of the block considered. Any truncated Taylor 
expansion, in this region, provides an increasingly poor approximation of tr p n A as 
the size of A increases. This is the region that we label 'non-perturbative'. 

The transition between these two regions turns out to be a smooth cross-over. 
The dependence of the threshold n c {L) on the size of the subsystem obeys a 
power law described in Eq. (|4.6p . We determined the critical exponents ruling 
this power law for the Ising model in both one and two dimensions. It would be 
interesting to investigate whether this cross-over is somehow related to the phase 
transition close to n=2 recently observed in [10] near d = 3. 

The above scenario has important consequences when considering the analyt- 
ical continuation involved in the computation of the entanglement entropy with 
the replica trick. This continuation is possible in one and two dimensions, as 
there is no phase transition between the perturbative region (where both ana- 
lytical results from QFT and numerical results from Monte Carlo algorithms are 
available) and the non-perturbative region (where the entanglement entropy is 
located) . However its form is not the naive expression obtained by promoting the 
simple expression obtained in the perturbative region for integer n to general a, 
as we have explicitly shown in Section [6] . 

From our analysis, the Renyi and Tsallis entropies seem to emerge as a privi- 
leged measures of entanglement in QFT. Although they lack the nice theoretical 
information properties of the entanglement entropy [28], their analytic proper- 
ties as functions of the a parameter seem to contain some new information on 
the entanglement. This has been considered recently in the context of quantum 
criticality [291 E01 ID] and topological order [51] . 

LT acknowledges the financial support from the QuSim group of the University 
of Queensland, the stimulating discussions with P. Corboz and the enlightening 
comments on entanglement measures of Prof. G. Vidal. 
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